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Abstract 

Background: Osteogenesis is a highly regulated developmental process and continues during the turnover and 
repair of mature bone. Runx2, the master regulator of osteoblastogenesis, directs a transcriptional program essential 
for bone formation through genetic and epigenetic mechanisms. While individual Runx2 gene targets have been 
identified, further insights into the broad spectrum of Runx2 functions required for osteogenesis are needed. 

Results: By performing genome-wide characterization of Runx2 binding at the three major stages of osteoblast 
differentiation - proliferation, matrix deposition and mineralization - we identify Runx2-dependent regulatory 
networks driving bone formation. Using chromatin immunoprecipitation followed by high-throughput sequencing 
over the course of these stages, we identify approximately 80,000 significantly enriched regions of Runx2 binding 
throughout the mouse genome. These binding events exhibit distinct patterns during osteogenesis, and are associated 
with proximal promoters and also non-promoter regions: upstream, introns, exons, transcription termination site regions, 
and intergenic regions. These peaks were partitioned into clusters that are associated with genes in complex biological 
processes that support bone formation. Using Affymetrix expression profiling of differentiating osteoblasts depleted of 
Runx2, we identify novel Runx2 targets including Ezh2, a critical epigenetic regulator; Crabp2, a retinoic acid signaling 
component; Adamts4 and Tnfrsf 1 9, two remodelers of the extracellular matrix. We demonstrate by luciferase assays that 
these novel biological targets are regulated by Runx2 occupancy at non-promoter regions. 

Conclusions: Our data establish that Runx2 interactions with chromatin across the genome reveal novel genes, pathways 
and transcriptional mechanisms that contribute to the regulation of osteoblastogenesis. 



Background 

The development of bone tissue, new bone formation in 
the adult, mineral homeostasis and maintenance of bone 
mass are mediated by cells of the osteoblast lineage. These 
bone forming cells progress through a highly regulated dif- 
ferentiation program, with each subpopulation of cells ac- 
quiring stage-specific phenotypes that are characterized by 
distinct profiles of expressed genes [1]. Osteoblast com- 
mitment and differentiation are dependent on the appro- 
priate expression of Runx2 (Runt-related transcription 
factor 2), the master regulator of bone formation [1-3]. 
Ablation of Runx2 in mice results in the absence of a 
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mineralized skeleton, and disruption of Runx2 function 
causes bone defects in the human disorder cleidocranial 
dysplasia [4,5]. Runx2 controls a complex gene-regulatory 
network during osteoblastogenesis [5,6]. It upregulates a 
variety of osteoblast lineage-specific genes, including Osx 
(osterix), Ocn (osteocalcin), and Bsp (bone sialoprotein), 
and represses the expression of non-osteoblast genes such 
as PPAR-y (peroxisome proliferator- activated receptor 
gamma) and MyoD (myogenic differentiation), which are 
required for adipogenic and myogenic commitment, re- 
spectively [7-9]. Studies have demonstrated that Runx2 
controls gene expression by interacting with multiple clas- 
ses of co-regulatory factors [1,3,10]. In addition to the 
traditional transcriptional mechanisms, several epigenetic 
mechanisms have been identified for Runx2-mediated 
gene expression. For example, Runx2 is retained on 
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mitotic chromosomes as an epigenetic bookmarking factor 
to maintain cellular identity after cell division [11]. A sig- 
nificant body of evidence substantiates the contribution of 
Runx2 towards the epigenetic regulation of gene expres- 
sion during osteoblast differentiation, through interactions 
with histone deacetylases [12], histone acetyltransf erases 
[13], and SWI/SNF complex components [14]. 

Runx2 has been shown to activate or repress gene expres- 
sion through binding to ds-regulatory DNA elements, the 
Runx motif (TGTGGT) located in or near gene promoter 
regions [15]. Genome-wide binding profiles of the tran- 
scriptional factor CTCF and lineage-specific transcription 
factors, such as PPAR-y, MyoD, and GATA3 (GATA bind- 
ing protein 3), via ChlP-Seq (chromatin immunoprecipita- 
tion followed by high-throughput sequencing) have helped 
to delineate ds-regulatory networks that are critical for cell 
lineage control in adipocytes, myocytes and T cells, respect- 
ively [16-19]. These studies have highlighted the regulatory 
importance of long-range interactions and binding of 
phenotypic transcription factors to non-promoter genomic 
elements. These findings have greatly expanded existing 
paradigms of transcriptional regulation. In contrast, the en- 
tire scope of Runx2 binding elements remains largely un- 
known, hindering a comprehensive understanding of the 
ds-regulatory network through which Runx2 regulates the 
transcription program for bone formation. 

Hence, we have characterized the genome-wide occu- 
pancy of Runx2 by ChlP-Seq in MC3T3-E1 preosteoblasts, 
a well-studied in vitro model for osteoblastogenesis [20]. 
Runx2 occupancy was determined at three hallmark stages 
of osteoblast differentiation: proliferation, matrix depos- 
ition, and mineralization. Analyses of these data demon- 
strated that Runx2 occupancy on genes is differentiation 
stage-dependent in both binding intensities and binding 
regions, indicating a shift in the regulatory mechanisms 
required for the entire program of osteoblastogenesis. By 
coupling genome-wide Runx2 binding with gene expres- 
sion profiling, we have identified new Runx2 targets that 
were validated for functional activities of Runx2 binding in 
both promoter and non-promoter regions. Our study of 
Runx2 genome-wide occupancy establishes a foundation 
for future investigation of the Runx2 -controlled regulatory 
network during bone formation and homeostasis. 

Results 

Runx2 dynamically occupies a wide range of genomic loci 
during osteoblastogenesis 

Runx2 is a known master activator of bone formation, 
but thus far only a small number of osteoblast-specific 
target genes have been characterized [2,5,6]. To identify 
genome-wide occupancy of Runx2 in osteoblast lineage 
cells, we performed ChlP-Seq using a Runx2- specific 
antibody at stages during the in vitro differentiation of 
MC3T3-E1 preosteoblasts (Figure 1). This in vitro cell 
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Figure 1 The differentiation stages of murine MC3T3-E1 

preosteoblasts used for profiling studies. (A) Staining for alkaline 

phosphatase (ALP) activity (upper panel) and mineralization (Von Kossa, 

lower panel) in MC3T3-E1 cells during three stages of differentiation: 

proliferation (day 0), matrix deposition (days 9 to 21), and mineralization 

(day 28). (B) Protein (left panel) and mRNA levels (right panel) of Runx2 

during osteoblastic differentiation of MC3T3-E1 cells. Histone H3 (H3) was 

used as loading control for western blotting. (C) Expression profile of 

osteoblast-related markers 0sx/Sp7, Collol (left panel), Akp2/Alpl (alkaline 

phosphatase liver/bone/kidney), Bsp/lbsp, and 0cn/Bglap2 (right panel) in 

MC3T3-E1 cells during differentiation. Relative mRNA levels (versus day 0) 

were determined by quantitative RT-PCR (reverse transcription PCR), 

normalized by Hprt 1 mRNA levels and plotted as mean values ± SEM 

(standard error of mean) from three independent biological replicates. 

The expression levels of the genes in (B,C) at days 9, 21, and 28 were 

significantly upregulated when compared to those at day 0 {P< 0.05, t-test). 
v J 

model recapitulates in vivo osteoblast differentiation and 
therefore was used for ChlP-Seq studies [21]. Alkaline 
phosphatase activity, an early osteoblast marker, in- 
creased as cells proceeded from proliferation (day 0) to 
matrix deposition (days 9 and 21) and decreased upon 
mineralization (day 28), visualized by Von Kossa staining 
(Figure 1A). Runx2 mRNA and protein levels signifi- 
cantly increased during the initial stage of differenti- 
ation. While the mRNA levels reached steady state, 
Runx2 protein levels declined during late mineralization 
(Figure IB). Osteoblast phenotypic markers, including 
the transcription factor Osx/Sp7 (osterix/Sp7 transcrip- 
tion factor), a marker of committed osteoprogenitors, 
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the extracellular matrix protein Collal (collagen type I 
alpha 1), and the specialized mineral binding proteins 
Bsp/Ibsp (bone sialoprotein/integrin-binding sialoprotein) 
and Ocn/Bglap2 (osteocalcin/bone gamma-carboxyglutamate 
(gla) protein 2), displayed expression patterns consistent with 
the progression of osteoblastogenesis [2] (Figure 1C). 

We next isolated and sequenced DNA from chromatin 
bound by Runx2. Sequence reads from Runx2 ChlPs and in- 
put controls were mapped to the mouse genome. Statisti- 
cally significant enrichments of Runx2 were identified by 
MACS (Model-Based Analysis of ChlP-Seq) [22] (Additional 
file 1). RefSeq annotations [23] were used to assign Runx2 
enrichments to categories of genomic locations (with non- 
overlapping definitions for transcriptional start site (TSS), 
promoter, exonic, intronic, transcription termination site 
(TTS), upstream and intergenic regions (see Figure 2 for 



details; Additional file 2). As Runx2 protein levels changed, 
the total number of Runx2 peaks changed accordingly 
(Figure 2 A, left panel). The overall distribution of Runx2 
binding among the categories of genomic locations was rela- 
tively unchanged during differentiation (Figure 2A, right 
panel). Among genomic locations, Runx2 occupancy at pro- 
moters showed the greatest variation from 17.6% at day 0 to 
8.8% at day 9 (Figure 2A, right panel). The majority of 
Runx2 binding occurred at intergenic and intronic regions 
(Figure 2A,B). However, when compared with a randomly 
sampled background distribution of genomic intervals 
(Additional file 3), Runx2 binding displayed preferential en- 
richment in a genie context, particularly at promoters and 
exons. One exception was under-represented Runx2 bind- 
ing at intergenic regions when compared to nonspecific or 
random binding (Figure 2B). This relative enrichment of 




Distance(bp) 



Figure 2 Genome-wide profile of Runx2 occupancy. (A) Distribution of Runx2 binding peaks across the mouse genome were classified into 

six categories of genomic locations: exon, intron, promoter (-1 kb to +150 bp of TSS), upstream (-1 kb to -20 kb from TSS), TTS region (-150 bp 

to -1 kb of transcription termination site), and intergenic region. Peak distribution was plotted at each time point by peak number (left panel) 

and by the percentage of total peaks (right panel). (B) Differential Runx2 enrichment in the six categories of genomic locations compared to the 

predicted Runx2 motif (see below), to random binding, or to binding of the transcription factor CTCF. Inset provides a magnified view of promoter 

and TTS regions. (C) The 500 most significant Runx2 peaks, based on MACS [22] significance (P< 1 x 10~ 10 ) were used for de novo motif discovery. A Runx2 

motif (position weight matrix, top) with strong statistical confidence (P=4.7 x 10~ 20 °) was determined using MEME (MEME suite version 4.7.0) [24]. The known 

Runx motif (MA0002.2, bottom) in JASPAR [26] was used for comparison using TOMTOM [24]. As shown in (B), the distribution of de novo Runx2 motifs 

among categories of genomic locations was determined using FIMO [24] at a significance threshold of P< 1CT 4 . (D) Probability plot of the distribution 

of Runx2 peaks indicating the distances of Runx2 motifs to the peak centers in the top, middle, and bottom third of Runx2 peaks (ranked by MACS 

scores) versus the probabilities of finding de novo Runx2 motifs at given positions relative to peak center. 
\ J 
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Runx2 occupancy in genie contexts was also observed in 
comparison with Runx2 motifs (Figure 2B). A de novo 
Runx2 motif (Figure 2C, top) was discovered using MEME 
[24] on the highest-ranked (by P-value) 500 ChlP-Seq 
peaks and was then scanned using FIMO [24] (version 
47.0) over the mouse genome. Notably, analysis by MEME 
did not discover a significantly enriched secondary motif 
associated with Runx2 peaks. Despite the nearly random 
distribution of Runx2 motifs throughout the genome, 
Runx2 occupancy in differentiating osteoblasts is character- 
ized by associations with promoters, exons, introns and 
other genie elements. These associations are perhaps due 
to epigenetic factors, including chromatin conformation 
and accessibility, along with co-factors, and suggest that 
the presence of a Runx motif does not necessarily indicate 
the physical association of Runx2. Interestingly, the distri- 
bution of Runx2 occupancy among classes of genomic ele- 
ments is similar to that of CTCF (Figure 2B), a ubiquitous 
transcription factor that exhibits a broad spectrum of DNA 
binding in many cell lines [25]. 

The de novo Runx2 motif (Figure 2C, top) was scanned 
over all ChlP-Seq peaks called by MACS [22] from ChlP- 
Seq reads collected from day 9 MC3T3 cells, and a histo- 
gram of the distance between the peak summit and the 
highest-scoring motif instance was collected. The distribu- 
tion in Figure 2D is characterized by a sharp peak for 
Runx2 motifs at the summits of ChlP-Seq peaks. Using 
TOMTOM [24] it was determined that the de novo Runx2 
motif was similar (E-value = 2.6 x 10" 4 , q-value = 2.6 x 10" 4 ) 
to the JASPAR MA0002.2 RUNX motif [26] (Figure 2C, 
bottom) and both motifs share the core TGTGGT se- 
quence with the known Runx2 binding consensus. 

Genome-wide Runx2 occupancy reveals distinct positional 
and temporal binding patterns 

We grouped Runx2 peaks based on the presence or ab- 
sence of Runx2 binding at specific time points during 
differentiation (Figure 3 A, off/on). The resulting seven 
distinct clusters reflected the dynamics of Runx2 bind- 
ing in relation to the progression of osteoblastogenesis 
(Figure 3A). In Figure 3B, the mean ChlP-Seq read dens- 
ities (that is, peak intensities) of the clustered Runx2 peaks 
were plotted to compare their relative enrichments. 

The seven clusters differed in the numbers and inten- 
sities of Runx2 peaks (Figure 3A). The largest group was 
cluster 6, which exhibited the presence of peaks primarily 
at day 9. This finding is consistent with day 9 representing 
committed osteoblasts with the highest amount of cellular 
Runx2 protein (Figure IB) and therefore the greatest num- 
ber of Runx2 peaks (Figure 3B). The strongest Runx2 peak 
intensities were found in cluster 1 reflecting regions that 
were bound by Runx2 constitutively throughout osteoblast 
differentiation (Figure 3B). The weakest Runx2 peak inten- 
sities among the three time points occurred at the peaks 



in cluster 7, with marginal enrichment of Runx2 at day 0. 
It is noteworthy that the peaks in cluster 4 have the sec- 
ond highest peak intensities at both days 9 and 28 (matrix 
deposition and mineralization stages; Figure 3B). Cluster 
5, like clusters 1 and 4, exhibited the highest peak inten- 
sities on day 28, indicating their importance in maintain- 
ing osteoblast phenotype. 

Runx2 binding in each cluster was further examined for 
distribution preferences of peaks in different genomic re- 
gions, in contrast to the genome-wide distribution of ran- 
dom 100 bp DNA fragments (detailed in Additional file 3; 
Figure 3C). The random DNA fragments (grey bars in 
Figure 3C) are distributed mainly in intergenic, intronic, 
and upstream regions. When compared to the random con- 
trol, we observed that the distribution of Runx2-enriched 
peaks was biased towards gene regions (exons, introns, pro- 
moters, TTS regions, upstream regions). In contrast, Runx2 
binding in the intergenic regions is lower than the ran- 
dom control (Figure 3C). In promoters and exons, all 
clusters showed the highest enrichment over random 
binding, suggesting strong regulation by Runx2 at these 
genomic regions. 

To further explore the relationship between Runx2 
peaks and peak-associated genes in osteogenic differenti- 
ation, we performed functional annotations for the peaks 
in the seven clusters using GREAT (Genomic Regions En- 
richment of Annotations Tool; Figure 3D). Gene Ontology 
(GO) terms associated with the largest clusters are shown 
in Figure 3D. Runx2 binding in cluster 1 yielded GO terms 
of general biological processes such as protein folding and 
RNA metabolism (detailed in Additional file 4). Cluster 4 
peaks at days 9 and 28 were frequently related to the GO 
terms of osteoblast differentiation, bone developmental pro- 
cesses, and osteogenic signaling pathways. These terms often 
included differentiation-related and well-known Runx2 tar- 
get genes, such as Runx2, Bsp/lbsp, and Osx/Sp7 (Additional 
file 4). Similarly, cluster 6 (day 9) peaks often associated with 
bone formation and extracellular matrix organization. The 
annotation of other smaller clusters is shown in Additional 
file 5. For examples, cluster 3 peaks were associated with 
apoptosis, programmed cell death, and DNA damage; clus- 
ter 7 containing peaks found only in day 0 was linked with 
negative regulation of cell cycle control and the phenotypes 
of non-osseous mesenchyme-derived cells (Additional files 4 
and 5). These functional annotations associated with Runx2 
peaks are generally consistent with the progression of 
MC3T3-E1 differentiation, supporting a temporal transcrip- 
tion network programmed by Runx2. 

Runx2 binding patterns at osteogenic genes are 
predictive of potential Runx2 targets 

To discover previously unknown Runx2 target genes, we 
first determined the Runx2 binding patterns of well-known 
Runx2 target genes found in differentiation cluster 4; for 
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Figure 3 Characterization of the temporal patterns of Runx2 binding associated with osteoblastogenesis. (A) Enriched regions of Runx2 
binding (MACS peaks, P< 1 x 10" 10 ) were grouped into seven clusters based on the presence (red) or absence (blue) of peaks at temporal stages 
of osteoblastic differentiation (left panel). In the right panel, each line illustrates intensity (base-2 logarithm of the sum of read counts per 10 
million reads (mean read density)) of one peak as a range from strong occupancy (brown) to weak occupancy (violet). (B) Runx2 binding 
intensities of Runx2 peaks from each cluster at days 0, 9, and 28 of differentiation plotted as mean ± SEM. (C) Distribution patterns of the Runx2 
peaks from each of the seven clusters (from (A)) segregated into six categories of genomic locations and compared to an equivalent number of 
random 100 bp regions (gray bar). (D) Representative Gene Ontology term annotation of cluster 4 (days 9 & 28) by GREAT. Each term is annotated 
with the observed (Obs.) number of Runx2 peaks and corresponding number of genes out of the total number of genes. 



example, Runx2, Osx/Sp7, and Ocn/Bglap2 (Figure 4A; 
Additional file 6), and Bsp/Ibsp (Additional file 7). For these 
genes Runx2 binding was distributed in promoters, the gene 
body (introns/exons) and sites distal from the gene body. 
Furthermore, Runx2 enrichment increased in the loci of 
these genes during osteoblast differentiation at matrix and 
mineralization stages. We then examined the genes associ- 
ated with cluster 1 peaks (Additional file 8), and identified 
genes, including Ezh.2 (enchancer of zeste homolog 2) 
(Figure 4B), with Runx2 binding profiles that displayed a 
ubiquitous but increasing level of Runx2 occupancy. Ezh2 is 
a component of PRC2 that epigenetically regulates gene ex- 
pression by methylating histone H3 lysine 27 and was 
recently found to be involved in commitment of mesenchy- 
mal stem cells towards the osteoblast lineage [27]. The 5' 
proximal promoter region of Ezh2 is bound by Runx2 
(Figure 4B; Additional file 6) and shows an increase in reads 
(occupancy) during differentiation. When Ezh2 mRNA 



levels were measured during osteoblast differentiation, we 
found that the highest expression occurred in proliferating 
MC3T3-E1 cells, when the lowest amount of Runx2 binding 
was observed (Figure 4C). The striking decrease in Ezh2 
mRNA levels with increased Runx2 binding suggests that 
the transcription of Ezh2 is potentially regulated by Runx2. 

Genes regulated by Runx2 exhibit distinct Runx2 binding 
profiles 

To characterize the Runx2 binding profiles in genes that 
are transcriptionally regulated by Runx2, we performed 
gene expression profiling at day 9 in MC3T3-E1 cells 
treated with control (Scr) short hairpin RNA (shRNA) or 
a Runx2-specific shRNA (shRunx2) that knocks down 
Runx2 expression. Runx2 protein levels decreased by 80% 
in cells treated with shRunx2 (Additional file 9). This 
knockdown in turn inhibited expression of differentiation 
marker genes and osteoblastogenesis as demonstrated by 
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Figure 4 Runx2 binding at known osteogenic genes and a potential target Ezh2, during osteoblastogenesis. (A) Runx2 binding at three 
known osteogenic genes, Runx2, Osx, and Ocn, during three stages of osteogenesis: proliferation (day 0); matrix deposition (day 9); and mineralization 
(day 28). Gene annotation follows standard gene prediction display conventions used by the UCSC genome browser (exons, solid boxes; introns, solid 
lines; direction of gene transcription, arrows). Positions of Runx2 peaks called by MACS (green bars) and Runx2 consensus motif (JGTGGT; solid black 
bar) are also depicted. (B) Ezh2 locus bound by Runx2 with increasing Runx2 enrichment over temporal stages of osteogenesis. (C) Ezh2 expression 
during osteogenesis (relative to day 0 levels) were determined by quantitative RT-PCR and normalized by Hprtl mRNA levels. Relative mRNA levels are 
plotted as mean ± SEM from three independent biological replicates. 



decreased Alp staining (Additional file 9). We found 159 
genes whose expression was responsive to Runx2 knock- 
down, with |log 2 (I s hR U nx2/Iscr)| > 1.5, and a false discovery 
rate (FDR) <0.05, where I is the measured normalized 



probe-set intensity. These genes included 115 up- and 44 
down-regulated genes (Figure 5A). The 15 genes most re- 
sponsive to Runx2 knockdown (Table 1) included the 
well-defined Runx2 targets Ocn/Bglap2, Bsp/Ibsp, and 
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Figure 5 Correlation of Runx2 occupancy and Runx2-reponsive genes identifies novel targets. (A) Expression levels of genes responsive to 
Runx2 silencing in differentiated MC3T3-E1 cells. Values are plotted as the log2(expression level) from shRunx2-expressing cells (vertical axis) versus control 
shRNA expression (horizontal axis). Each point represents mean mRNA expression level from three independent biological replicates. Several representative 
genes are labeled. Diagonal lines demarcate the threshold for significant increase or decrease (>1 .5-fold) in expression. (B-D) Runx2 peaks associated with 
upregulated (Up), downregulated (Down), or unchanged (Non-responsive) gene expression upon Runx2 knockdown were compared by: average peak 
number per gene (B), and peak distribution (C) and fold change of peak signals (d9 versus dO) across genomic locations (D) with shRunx2 non-responsive 
genes as a control. In (B), three groups were compared to non-responsive genes: all peaks in shRunx2-regulated genes (All peaks), all peaks present at day 
9 in shRunx2-regulated genes (D9 peaks), and all peaks present at day 0 in shRunx2 regulated genes (DO peaks). Values are mean ± SEM (B 7 D) and statistical 
significance {*P< 0.01, **P< 0.05) determined by Mann-Whitney test (B 7 D) or Fisher's exact test (C). (E) Runx2 enrichment across gene bodies (±10 kb) of 
genes downregulated (Down), upregulated (Up), and unchanged (Non-responsive) by shRunx2 treatment at day 9. Mean signal ratios (IP/Input) at each 
genomic region were determined using PeaksToGenes. Error bars represent SEM. (F) Mean phyloP conservation scores of Runx2 motifs associated with 
genes significantly (fold change >1.5, FDR <0.05) downregulated (Down), upregulated (Up), or unchanged (Non-responsive) upon shRunx2-treatment. 
Conservation of Runx2 motifs was compared between shRunx2 downregulated and upregulated genes and length-matched non-responsive genes 
(Non-responsive, Down and Up, respectively) and statistical significance determined (**P<0.05) by Kolmogorov-Smirnov test. 
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Table 1 Top 15 genes upregulated or downregulated during 
shRunx2 treatment when compared to scramble shRNA 



Down 




Up 




Symbol 


Fold change 


Symbol 


Fold change 


Phex 


0.26 


l830012O16Rik 


4.73 


H19 


0.36 


lfi27l2o 


3.83 


Ibsp 


0.46 


Oos2 


3.78 


Mmpl3 


0.46 


D14Ertd668e 


3.62 


Mid2 


0.49 


Ifit3 


3.37 


Edil3 


0.50 


Gm12250 


3.27 


TnfrsfW 


0.53 


Ifitl 


3.25 


Bglop2 


0.56 


Ly6c1 


3.24 


C87414 


0.56 


Usp18 


3.09 


Gpr116 


0.56 


BC023105 


3.05 


Pdliml 


0.57 


ligpl 


3.05 


Hs3st3a1 


0.57 


Ddx60 


2.88 


Id3 


0.58 


Gm4951 


2.87 


Choc! 


0.58 


Ly6o 


2.87 


1133 


0.58 


Ly6c2 


2.85 



Mmpl3, which are known to be activated by Runx2 [5]. 
Notably, the genes most upregulated by shRunx2 have not 
been characterized as Runx2 targets (Table 1) except for 
Uspl8, which encodes a protein involved in the ubiquitin 
degradation pathway [28] . 

We characterized the relationship between Runx2 oc- 
cupancy and genes affected by shRunx2 knockdown, 
compared to non-responsive genes. The number of 
peaks, genomic distribution of peaks, and fold change of 
peak signals were compared among the gene groups at 
days 0 and 9 (Figure 5B-D). Genes that were downregu- 
lated by shRunx2 at both days 0 and 9 had more Runx2 
peaks when compared to the control non-responsive 
genes (Figure 5B). It should be noted, however, that 
shRunx2 downregulated genes are longer than the up- 
regulated genes (Additional file 10). This finding was 
also observed when we included the Runx2 binding pro- 
files at day 28 for Runx2-regulated genes (Additional file 
10). In contrast, the shRunx2 upregulated genes that ap- 
peared on day 9 had fewer Runx2 peaks compared to 
control genes (Figure 5B). There were also significant 
differences in overall peak distribution between shRunx2 
responsive genes and control (Figure 5C). Genes that were 
downregulated tended to exhibit more intronic and inter- 
genic enrichment of Runx2 peaks, while shRunx2 upregu- 
lated genes were strongly enriched in intergenic but 
reduced in intronic binding. We further examined the dis- 
tribution of Runx2 peaks in up- and down-regulated genes 
as a function of changes in Runx2 binding during differen- 
tiation from days 0 to 9 (Figure 5D). The fold change in 
day 9/day 0 peak signals showed increased Runx2 binding 



predominantly at upstream and promoter regions for the 
shRunx2 downregulated genes; in shRunx2 upregulated 
genes, Runx2 binding diminished in the introns and 
showed no significant change in promoter regions. There- 
fore, shRunx2 downregulated and upregulated genes ex- 
hibited distinct preferences for Runx2 binding in genomic 
loci as reflected by peak distributions in Figure 5B,C, and 
in relation to differentiation (Figure 5D). 

To complement the above analysis examining the 
genome- wide distribution of Runx2 responsive peaks, we 
used the PeaksToGenes program [29] to determine the 
enrichment of day 9 Runx2 signals at defined intervals 
within and surrounding the gene bodies (Figure 5E). For 
shRunx2 downregulated genes, Runx2 binding had the 
strongest enrichment surrounding TSSs, which includes 
proximal promoter, 5 ' UTR, exons and introns all within 
five contiguous deciles (Figure 5E; Additional file 11). 
This analysis demonstrated that, of genes affected by 
shRunx2, 66.1% of upregulated genes have a low level 
of Runx2 binding (cluster V in Figure S5B in Additional 
file 11) but most downregulated genes (69.1%) have 
higher level intensities of Runx2 binding (Figure S5A,B 
in Additional file 11). This finding shows that Runx2 re- 
sponsive genes at day 9 (shRunx2/Scr) are primarily reg- 
ulated by Runx2 surrounding TSSs. 

As another computational analysis, EMBER (Expect- 
ation Maximization of Binding and Expression pRofiles) 
[30] was used to relate measured changes in gene ex- 
pression to the spectrum of Runx2 occupancy observed 
during osteoblast differentiation (Figure 3B,C; Additional 
file 12). In analogy with discovering a sequence motif 
from a collection of functionally related DNA sequences, 
EMBER optimizes an 'expression pattern' from a collec- 
tion of genes related by patterns of transcription factor 
binding and uses this motif to determine which genes 
are regulatory targets of the transcription factor (details 
of EMBER are summarized in Additional file 3). Using 
this approach, we discovered and compared expression 
patterns from different sets of Runx2 binding regions 
(Figure 3). The Runx2 peaks were partitioned into 42 
subsets of Runx2 binding regions (7 clusters with 6 gen- 
omic location categories) and it was observed that all 
groups of Runx2 binding show a correlative relationship 
to gene expression during osteoblast differentiation 
(Additional file 12). We noted, however, that intronic 
Runx2 binding regions - for all time-dependent patterns 
of Runx2 occupancy - are less informative than other 
groups, indicating that intronic Runx2 binding may be 
less useful than binding at other class elements as a pre- 
dictor of gene regulation. 

Finally, downregulated and upregulated genes were 
compared on the basis of their evolutionary conservation 
(Figure 5F). Functional genomic elements are often char- 
acterized by conservation, which has been used to guide 
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the prediction of transcription factor binding sites 
[31,32]. This finding has been shown to distinguish func- 
tionally verified from un-verified binding sites in a large- 
scale study of transcription factor binding site function 
on human promoters [33]. For each ChlP-Seq peak, a 
Runx motif was used to identify the single most likely 
Runx2 binding site and the mean phyloP score [34] (for 
conservation among 30 vertebrate species) across the 
binding sites was computed. For each gene, phyloP 
scores were averaged among all ChlP-Seq peaks that 
were associated with individual genes to give an average 
measure of Runx2 conservation. We found that genes 
downregulated by shRunx2 had Runx2 binding se- 
quences that were more conserved than those in upreg- 
ulated or non-responsive genes (Figure 5F), suggesting 
that Runx2 regulation may be more evolutionarily con- 
served in genes that are activated by Runx2 during 
osteoblastogenesis. 

Taken together, the complementary methods described 
above suggest that Runx2 employs different mechanisms to 
regulate gene expression: 1) shRunx2 downregulated genes 
show increased Runx2 binding at promoter and far up- 
stream regions; and 2) based on EMBER, intronic binding 
does not imply Runx2-mediated gene regulation to the same 
degree as Runx2 binding at promoter, exon and upstream 
regions. Thus, during the normal course of osteoblast differ- 
entiation, Runx2-activated genes (shRunx2 downregulated) 
are regulated through both promoter and non-promoter re- 
gions; and regulation of Runx2 repressed genes (shRunx2 
upregulated) also occurs in promoter and other genomic re- 
gions, but with less Runx2 binding. 

Novel genes regulated by Runx2 through distinct 
regulatory elements 

Among potential Runx2 targets, we identified Tnfrsfl9 
(tumor necrosis factor receptor superfamily, member 
19), which is involved in bone formation as a Wnt- 
responsive regulator of mesenchymal stem cell com- 
mitment to osteoblastic lineage [35]. Tnfrsfl9 exhibited 
enrichment of Runx2 binding to the promoter as well as 
intronic regions (Additional files 6 and 13). During dif- 
ferentiation of osteoblasts, Tnfrsfl9 mRNA levels in- 
creased dramatically more than 100-fold from day 0 to 
day 9 (Figure S7A in Additional file 13). Consistent with 
Affymetrix data, depletion of Runx2 resulted in signifi- 
cant decreases in Tnfrsfl9 expression, indicating direct 
Runx2 regulation (Figure S7B in Additional file 13). We 
found that Runx2 occupancy was increased at intronic 
and promoter regions of the Tnsrsfl9 locus during osteo- 
blast differentiation (Figure S7C in Additional file 13). 

Adamts4 and Crabp2 were selected for further ana- 
lyses based on their responsiveness to shRunx2 (with 
fold change cutoff of 1.3; Additional file 14) and their 
Runx2 binding patterns predominantly in non-promoter 



regions. To test the functionality of Runx2 binding at these 
putative regulatory regions, we cloned non-promoter Runx2 
binding regions and measured transcriptional activity 
by luciferase reporter assay. Adamts4 is expressed in 
osteoblasts and osteocytes and encodes an enzyme that 
degrades aggrecan [36,37]. Runx2 exhibits multiple peaks 
across the Adamts4 locus, with increased occupancy 
during osteoblast differentiation (Figure 6A; Additional 
file 6). Adamts4 expression during osteoblast differenti- 
ation was increased at day 9 and remained steady to day 
28 (Figure 6B). Knockdown of Runx2 suppressed the ex- 
pression of Adamts4 (Figure 6C) by 40%, supporting 
Adamts4 as a direct target of Runx2 during osteoblasto- 
genesis. We characterized the functional activity of two 
prominent Runx2 binding regions, peak A in intron 1 and 
peak B in the last exon (Figure 6A). The peak A region 
increased luciferase activity in MC3T3-E1 cells by over 
20-fold. In contrast, peak B functioned as a suppressor 
of luciferase activity (Figure 6D). These results indicate 
that the two Runx2 sites can function as a positive and 
negative regulator of Adamts4; however, the weaker, 
negative regulation by peak B on the luciferase reporter 
may be due to the lack of a native chromatin context. 
The large increase in luciferase activity observed from 
peak A is consistent with a previous study that demon- 
strated Runx2 upregulates ADAMTS4 in human SW1353 
chondrosarcoma cells [38]. Here we established in osteo- 
blasts that Runx2-mediated upregulation of Adamts4 can 
occur at non-promoter regulatory elements, and is not re- 
stricted to the proximal promoter region as previously 
shown [38]. 

Crabp2 is a cytoplasmic retinoic acid binding protein 
previously reported to be upregulated during osteoblasto- 
genesis [39] . We observed that Runx2 constitutively occu- 
pies the Crabp2 locus in the first intron, while binding 
increases upstream and downstream of the Crabp2 gene 
body during differentiation (Figure 6E; Additional file 6). 
Crabp2 was upregulated during MC3T3-E1 cell differenti- 
ation and knockdown of Runx2 decreased Crabp2 mRNA 
level (Figure 6F). In proliferating cells, peak regions C, D 
and E demonstrated minimal luciferase reporter activ- 
ity. However, in differentiating cells, all peaks exhibited a 
significant increase of luciferase activity, with peak region D 
showing a more than 30-fold activation (Figure 6G). The 
knockdown of Runx2 reduced luciferase activity (Figure 6H), 
further supporting the function of these regions in mediat- 
ing Runx2 regulation of Crabp2. 

These findings of novel genes bound and regulated by 
Runx2 through different types of genomic elements sup- 
port an emerging concept that non-promoter elements 
can regulate gene transcription. Our results also indicate 
that Runx2 mediates complex fine-tuning of gene expres- 
sion in osteoblasts by both activating and repressing regu- 
latory elements that are located in non-promoter regions. 
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(See figure on previous page.) 

Figure 6 Non-promoter association of Runx2 regulates novel targets Adamts4 and Crabp2. (A) Increasing Runx2 enrichment was observed 
in the first intron (peak A, boxed region) and last exon (peak B, boxed region) of the Adomts4 locus during osteogenic differentiation. (B) Adomts4 
mRNA levels (normalized to Hprtl) were significantly upregulated (P< 0.05) in differentiating MC3T3-E1 cells (days 9 to 28) when compared to proliferating 
cells (day 0). (C) Runx2 knockdown (shRunx2) decreases Adamts4 expression, compared to a scrambled shRNA (Scr) control {**P < 0.05). (D) DNA sequences 
identical to peak regions A and B were cloned into individual pGL2-SV40-Luc reporters and relative luciferase activity was measured in transfected MC3T3-E1 cells 
and significant increases and decreases {*P< 0.01) in luciferase activity were observed for peak A and B reporters, respectively. (E) Increasing Runx2 enrichment 
was observed up- (peak C) and downstream (peaks D and E) of the Crabp2 locus during osteogenic differentiation. (F) Crabp2 expression increases during 
differentiation. Knockdown of Runx2 (by shRunx2) significantly reduces (**P<0.05) the endogenous expression of Crobp2 (right panel). (G) DNA sequences 
identical to peak regions C, D and E were cloned into individual pGL2-SV40-Luc reporters and relative luciferase activity was measured in transfected MC3T3-E1 
cells at days 0 and 7 {*P < 0.01). (H) Runx2 knockdown by inducible shRNA results in a significant decrease {**P < 0.05) of luciferase activity mediated by peak C 
and E regions and a downward trend (P = 0.08) in luciferase activity regulated by peak D region. Statistical significance for all experiments was determined by 
Student's t-test from mean ± SEM from three biological replicates. 



Discussion 

Through comprehensive genomic analysis of Runx2 by 
ChlP-Seq, we describe widespread Runx2 binding through- 
out the genome of differentiating osteoblasts. In addition 
to Runx2 interaction at promoters, we find Runx2 binding 
in non-promoter regions regulating novel targets that are 
silenced or expressed at different stages of osteoblast dif- 
ferentiation. Runx2 peaks cluster into temporal and func- 
tional categories associated with genes in a broad range of 
cellular programs, including bone development, negative 
regulation of proliferation, active matrix formation and 
mechanisms for mineral deposition that reflect the progres- 
sion of osteogenesis. Our data have identified new Runx2- 
regulated genes (Tnfrsfl9, Adamts4, Crabp2, and Ezh2) that 
have established roles in bone formation, and more import- 
antly, have extended the understanding of Runx2-mediated 
gene regulation to a broader range of cellular functions dur- 
ing osteoblast differentiation. 

Runx2 binding patterns identify stage-dependen 
t osteogenic programs 

Time-dependent Runx2 binding patterns underlie the dy- 
namic gene regulation by Runx2 during osteohistogene- 
sis. We identified a large number of peaks, consistent with 
the increasing protein levels of Runx2 from the early 
osteoprogenitor to the mature osteoblast/osteocyte. The 
results from binary clustering, together with subsequent 
GO term analyses by GREAT, identified many categories 
associated with skeletal development and bone homeosta- 
sis. These findings support our initial hypothesis that dis- 
tinct binding patterns of Runx2 at different stages of 
osteoblastogenesis have novel functional implications. 

Our clustering analysis partitioned Runx2 peaks into two 
main categories: less variable steady-state binding (cluster 
1, ubiquitous peaks) and more dynamic binding groups 
(clusters 2 to 7, stage-specific clusters). Steady-state binding 
of Runx2 persists from proliferative pre-osteoblasts to dif- 
ferentiating osteoblasts with binding signals plateauing dur- 
ing matrix deposition and mineralization stages. Peaks in 
this category represent genes related to housekeeping 
processes such as protein folding, negative regulation of 



mitotic cell cycle, and mRNA catabolism and processes not 
previously related to Runx2 (Additional file 4). The GO 
terms associated with stage-specific clusters include nega- 
tive regulators of other cell lineages (that is, fat and smooth 
muscle cells) as well as positive regulators of osteogenesis. 
Thus, dynamic Runx2 binding primes, enhances and stabi- 
lizes the osteoblast phenotype as well as suppresses non- 
osteoblast lineages. Many Runx2 bound genes in cluster 4 
(days 9 & 28, differentiation) have been demonstrated to 
contribute to in vivo bone formation [40]. Thus, the gen- 
omic profiling of Runx2 binding in our in vitro model sys- 
tem is consistent with the known properties of Runx2 
in bone formation. More importantly, our profiling re- 
veals pathways previously unknown to be controlled by 
Runx2 underlying biological mechanisms of general cellular 
processes. 

Runx2 binding functions at both non-promoter and 
promoter regions 

Only a small proportion of sequence-specific transcription 
factors, such as Myc, have narrow distributions of binding 
around proximal promoters of genes [41-44]. In contrast, 
non-promoter binding is recognized for other transcrip- 
tion factors, including STAT1, RUNX1, ERa, CTCF, and 
HNF4a [41,45-49]. In a previous study, overexpression of 
Runx2 in prostate cancer cells revealed extensive non- 
promoter binding [50]. In our study, endogenous Runx2 
binding across the genome was characterized during os- 
teoblastogenesis. We found that over 70% of Runx2 occu- 
pancy was localized to non-promoter regions (intergenic, 
intron, exon, TTS, and upstream regions that consti- 
tute the bulk of the genome) during the differentiation of 
osteoblasts. From days 0 to 9 there was a two-fold in- 
crease in the number of non-promoter peaks, indicating a 
functional association with the differentiation process. 
Runx2-dependent regulation through non-promoter peaks 
around the Adamts4 and Crabp2 genes provided direct 
evidence that non-promoter binding of Runx2 controls 
gene expression. 

Although we have demonstrated the importance of 
non-promoter Runx2 binding events, Runx2 peaks at 
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promoter regions have critical regulatory roles as well. 
In our data, Runx2 occupancy has the highest enrich- 
ment at promoter regions when compared with other 
genomic locations. Examples of this regulatory mode 
can be seen in some well-characterized Runx2 targets 
such as Bsp and Ocn, in line with established evidence 
that Runx2 can regulate these genes in a promoter- 
dependent manner [51-54]. The genes downregulated 
upon Runx2 silencing also displayed a clear enrichment 
of Runx2 signal at promoter regions (Figure 5B), exem- 
plified in Tnfrsfl9. 

Gene expression in different biological settings is influ- 
enced by higher order three-dimensional chromatin 
complexes that involve looping of promoter and non- 
promoter elements, blurring the distinction of defined 
regulatory regions [55,56]. It is plausible that some 
Runx2 peaks in promoter and non-promoter regions 
may serve as nucleation sites for modifications of chro- 
matin structures necessary for gene expression. Further- 
more, transcription factors can interact with RNA 
polymerase II, CTCF, and other factors via higher-order 
chromatin conformations [55,57,58]. During osteoblasto- 
genesis, Runx2 exhibited a distribution pattern among 
genomic elements similar to that of CTCF, suggesting 
that, like CTCF, Runx2 may have a functional role that 
extends beyond direct regulation of transcription. Con- 
sistent with this idea, Runx2 is able to form discernible 
foci associated with the nuclear matrix: a nuclear frame- 
work for organizing higher order chromatin structures. 
Runx2 truncated of the NMTS (nuclear matrix-targeting 
signal) domain results in diminished nuclear matrix as- 
sociation and disrupted expression of Runx2 target genes 
[59]. It is also noteworthy that the genes upregulated 
upon Runx2 knockdown have preferential Runx2 bind- 
ing in intergenic regions, indicating that distal elements 
may have a regulatory role through long-range interac- 
tions. In addition, Runx2 binding at distal regions may 
contribute to chromatin remodeling through interacting 
with chromatin-modifying enzymes, as has been well 
documented at regulatory elements in osteoblasts and 
other cell models [12,13]. 

The complexities of Runx2 binding and transcriptional 
regulation 

Runx2 displays complex binding patterns similar to other 
lineage-specific transcription factors, such as PPAR-y, 
MyoD, and GATA3 [16-19]. By systematic annotation of 
Runx2 peaks, multiple integrative analyses of gene expres- 
sion combined with Runx2 binding profiles and direct 
experimental validation of individual targets, we defined 
Runx2 binding with biological outcomes during osteoblast 
differentiation. These analyses revealed that, for a small 
set of genes, the enrichment and binding patterns of 
Runx2 were indicators of gene expression. Genes that have 



decreased expression in the absence of Runx2 (by 
shRunx2 treatment) have a greater number of Runx2 
peaks and greater fold change of peak signal at promoter 
and non-promoter regions when compared to non- 
responsive and gene length-matched controls. In contrast, 
genes upregulated by Runx2 knockdown tend to have 
fewer Runx2 peaks and smaller relative fold changes in 
peak signals when compared to controls. One notable 
finding that arose from our analysis was that genes that 
were downregulated in the absence of Runx2 had both a 
greater evolutionary conservation of Runx2 binding sites 
and tended to be longer than non-responsive and upregu- 
lated genes. It is unclear why this particular subset of 
genes (that is, genes normally activated by Runx2) would 
retain these features throughout evolution; however, this is 
an interesting point for future investigations. Although we 
demonstrated that Runx2 binding influences gene expres- 
sion, a proportion of Runx2 peaks were found to have no 
direct function in transcriptional control of genes. This may 
be due to efficient but not complete knockdown of Runx2 
by viral-mediated shRNA. Alternatively, similar to other 
transcription factors, many binding regions were found to 
be nonfunctional in transactivating luciferase reporters [33]. 
This finding suggests that binding of transcription factors 
may have a distinct function other than direct control of 
gene expression, consistent with previously described non- 
transcriptional functions and genome-organizing capabilities 
ofRunx2 [59,60]. 

Conclusions 

Our findings provide a new level of understanding of the 
Runx2-mediated transcription program as defined by 
genome-wide Runx2 binding essential for osteoblasto- 
genesis. Our data support that Runx2 functions at pro- 
moter and non-promoter regions at both previously 
known and novel targets. The impact of our study exam- 
ining the global occupancy of endogenous Runx2 in dif- 
ferentiating osteoblasts sets a framework for novel 
mechanisms underlying bone biology. 

Materials and methods 

Cell culture 

The calvaria-derived preosteoblast cell line MC3T3-E1 
(Subclone 4) was obtained from ATCC (Manassas, VA, 
USA) and maintained in ascorbic acid-free alpha-MEM 
(Hyclone, Novato, CA, USA) supplemented with 10% fetal 
bovine serum (FBS; Hyclone), 2 mM L-glutamine, 100 U/ml 
of penicillin and 100 (ig/ml streptomycin (Pen/Strep, Invi- 
trogen, Carlsbad, CA,). To induce osteogenic differentiation, 
complete alpha-MEM was supplemented with 280 uM 
ascorbic acid and 10 mM beta-glycerophosphate (Sigma 
Aldrich, St. Louis, MO, USA). Cells were maintained at 37°C 
in a humidified 5% C0 2 environment and media replaced 
every 2 to 3 days for the duration of all experiments. 
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RNA isolation and quantitative PCR 

Total RNA was isolated using Trizol reagent (Invitrogen) 
according to the manufacturers' specifications. Total cel- 
lular RNA treated by DNasel (Zymo, Irvine, CA, USA) 
was primed with random hexamers and reverse tran- 
scribed into cDNA using Superscript First-strand cDNA 
Synthesis Kit (Invitrogen) according to the manufac- 
turers instructions. Gene expression was determined by 
quantitative real time PCR (qPCR) using iQ™ SYBR Green 
PCR Master Mix (BioRad, Hercules, CA, USA) in an ABI 
Prism 7300 thermocycler (Applied Biosystems, Foster 
City, CA, USA). For each gene, the expression level was 
normalized to that of Hprtl using 2" AACt method. Experi- 
ments were performed in triplicate and results are pre- 
sented as mean values ± SEM. Primers for qPCR reactions 
were designed by FoxPrimer [61,62], and are available in 
Additional file 15. 

Runx2 knockdown and gene expression profiling 

Lentiviruses carrying Runx2-shRNA and previously de- 
scribed control Scramble-shRNA [63] were used to in- 
fect MC3T3-E1 cells. Infected cells were subsequently 
detected by green fluorescent protein, sorted and grown 
to 90% confluency followed by osteogenic differentiation 
for 9 days. Knockdown experiments were performed in 
three biological replicates. 

Microarray samples were handled following the manu- 
facturers' recommended protocols (Affymetrix, Santa 
Clara, CA, USA). Briefly, RNA isolated from MC3T3-E1 
cells (day 0 and 9 scramble shRNA, and day 9 Runx2 
shRNA) were reversely transcribed into cDNA using WT 
Expression Kit (Ambion, Austin, TX, USA), labeled and 
fragmented with GeneChip WT Terminal Labeling and 
Controls Kit. Labeled cDNAs were then hybridized to 
GeneChip Mouse Gene 1.0 ST Array rev.4 using a Gene- 
Chip Hybridization, Wash, and Stain Kit. Hybridization 
signals were obtained by GeneChip Scanner (Affymetrix). 
Microarray data were analyzed using Bioconductor (ver- 
sion 2.11) packages affy and limma in R (version 2.15.1) 
[64-67]. Briefly, after performing RMA (robust multichip 
average) normalization of microarray expression levels 
and filtering, differential expression was detected using a 
Bayesian moderated £-test. The Benjamini-Hochberg FDR 
[68] was applied to correct for multiple testing. The Affy- 
metrix expression profiles were annotated to RefSeq genes 
[23]. These expression data have been deposited in the 
Gene Expresion Omnibus (GEO) database under acces- 
sion number GSE53982. 

Western blot 

Nuclei extracts were prepared from MC3T3-E1 cells 
using a protocol modified from Dignam et al [69]. Pri- 
mary antibodies and dilutions were: mouse anti-Runx2 
monoclonal (Clone 8G5, MBL International, Woburn, 



MA, USA; 1:1,000); rabbit anti-H3 3H1 monoclonal 
(Cell Signaling, Danvers, MA, USA; 1:2,000). HRP- 
conjugated secondary antibodies and dilutions were: goat 
anti-mouse IgG (Santa Cruz, Dallas, TX, USA; 1:3,000); 
goat anti-rabbit IgG (Santa Cruz; 1:3,000). Detection of 
HRP was performed using a Western Lightning Plus Kit 
(Perkin Elmer, Waltham, MA, USA) followed by expos- 
ure onto Biomax Light film (Kodak, Rochester, NY, 
USA). 

Chromatin immunoprecipitation and high-throughput 
sequencing 

At days 0, 9, and 28 of differentiation, approximately 1 x 
10 8 MC3T3-E1 cells were washed with PBS (phosphate- 
buffered saline) and then fixed on a plate with 1% for- 
maldehyde for 8 minutes to crosslink DNA-protein com- 
plexes. The fixed cells were washed with ice-cold PBS, 
harvested, and pelleted. Nuclei extraction was performed 
using a protocol modified from Dignam et al [69]. Iso- 
lated nuclei were sonicated using a Misonix S-4000 
ultrasonic sonicator to obtain sheared chromatin ranging 
from 0.2 kb to 0.6 kb. Sheared chromatin was used 
for immunoprecipitation with Runx2 antibody (M-70, 
Santa Cruz) [70] or immunoglobulin G (IgG) (12-370, 
Millipore, Billerica, MA, USA) followed by purification 
using Protein-G Dynabeads (Invitrogen). Precipitated chro- 
matin was washed with solutions of increasing salt con- 
centration, and eluted and subsequently uncrosslinked at 
65°C. DNA was recovered by phenol/chloroform/isoa- 
myl alcohol extraction followed by ethanol precipita- 
tion. Libraries of purified DNA were generated using 
Illumina SR adapters (Illumina, San Diego, CA, USA) 
following the manufacturers manual. DNA libraries 
were selected for inserted fragments of 200 ± 50 bp, 
and single-end 36 base reads were generated on an Illu- 
mina Genome Analyzer II at the UMASS Deep-Seq 
core facility. Base calls and sequence reads were gener- 
ated by Illumina CASAVA software (version 1.6; Illu- 
mina). Two independent biological repeats of Runx2 
ChlP-Seq libraries were prepared for each time point, 
and two Input libraries were prepared with sonicated 
DNA from day 9 MC3T3-E1 cells. 

Analysis of ChlP-Seq data 

Single-end 36 base sequences from Runx2 ChlP-Seq and 
input libraries were mapped to the mouse genome (as- 
sembly mm9) using Bowtie (version 0.12.8) [71]. Runx2 
peaks and read counts were determined by MACS (ver- 
sion 1.4.1) [22] using default settings. Runx2 peaks that 
were significant at the P < 10" 10 level were retained for 
subsequent analysis. For each of the three time points in 
our data, read counts were normalized to 10 million 
reads. The UCSC genome browser [72] was used to 
visualize Runx2 peaks. 
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Runx2 binding regions were classified on the basis of 
genomic location categories and annotated to known 
RefSeq genes [23] . Runx2 peaks were grouped into seven 
clusters based on the presence or absence of a peak. 
Runx2 peaks in each cluster were analyzed for GO terms 
using GREAT (version 2.0.2), using default association 
rules between ChlP-Seq peaks and annotated genes [73]. 
The conservation of Runx2 motifs associated with the 
genes responsive to shRunx2 were determined using 
phyloP [34] . Statistical significance was determined using 
Kolmogorov-Smirnov test. 

PeaksToGenes [29,62] (Additional file 3) was used to 
test Runx2 binding in relation to the genes from micro- 
array analysis. PeaksToGenes defines genomic intervals 
relative to all RefSeq genes, and in each window uses a 
non-parametric Wilcoxon rank sum test to calculate the 
probability and binding frequency. The comparisons 
were individually made between each responsive group 
and the non-responsive group to Runx2 shRNA. Runx2 
binding and Runx2-mediated transcription control were 
also evaluated by EMBER [30]. Analogous to discovering 
a sequence motif from functionally related DNA se- 
quences [24], EMBER optimizes an expression pattern 
from a collection of genes' expression data related by 
profiles of transcription factor binding and uses this in- 
formation to determine which genes are potential regu- 
latory targets of the transcription factor. For a detailed 
description of the PeaksToGenes and EMBER analysis, 
please refer to Additional file 3. 

Functional genomic elements can be characterized by 
evolutionary conservation and have been shown to dis- 
tinguish functionally verified from unverified transcrip- 
tion factor binding sites on human promoters [33]. For 
each ChlP-Seq peak (combined from day 0, 9, and 28 
datasets), a Runx motif was used to identify the single 
most likely Runx2 binding site and the mean phyloP 
score [34] (for conservation among 30 vertebrate spe- 
cies) across the binding sites was computed. In order 
to compare conservation among genes that fell into 
different regulatory groups (shRunx2 downregulated, up- 
regulated and non-responsive genes) based on our ex- 
pression microarray measurements, phyloP scores were 
averaged among all ChlP-Seq peaks that were associated 
with individual genes to give an average measure of 
Runx2 conservation for each gene. To minimize the 
effect of gene length on Runx2 binding analysis, we 
compared the conservation of shRunx2-downregulated and - 
upregulated genes to randomly sampled, length-matched 
non-responsive genes. Statistical significance was deter- 
mined by pair-wise conservation comparisons using 
Kolmogorov-Smirnov test. 

The raw sequences and peak-related files in BED and WIG 
formats representing processed data have been deposited in 
the GEO database under accession number GSE54013. 



Luciferase assays and plasmid reporters 

Selected Runx2 peak regions were cloned with Mlul 
combined with Bglll or Xhol into a pGL2-SV40-Luc re- 
porter (Promega, Fitchburg, WI, USA). Primers used in 
cloning are listed in Additional file 16. A pGL2-SV40- 
Luc reporter with minimum SV40 promoter was used as 
mock control. Reporter plasmids with Runx2 peak regions 
and pGL2-SV40-Luc empty vector were co-transfected 
with pcDNA3.1-Runx2-WT or pcDNA3.1-EV into MC3T3 
cells. Transfected cells were then differentiated for 7 days 
before luciferase activities were determined by Dual-Glo 
Luciferase Assay Kit (Promega). A second set of luciferase 
assays was done in differentiating MC3T3-E1 cells stably 
infected with a doxycycline-inducible pLKO-puro-Tet-on- 
Rx2shRNA lentiviral vector. This vector was constructed 
by re-cloning a previously validated Runx2 shRNA se- 
quence [63] to Tet-pLKO-puro plasmid (catalog number 
21915, Addgene, Cambridge, MA, USA). Luciferase activ- 
ities from MC3T3-E1 cells treated with or without 2.5 
ug/ml doxycycline were examined at day 7 after differenti- 
ation (Figure 6). 

Additional files 



Additional file 1: Table S3. Summary of MC3T3 Runx2 ChlP-Seq. This table 
lists the read numbers and genome coverage of Runx2 ChlP-Seq libraries. 

Additional file 2: Table S4. Distribution patterns of Runx2 peaks across 
genomic locations. This table is related to Figure 2. 

Additional file 3: Detailed description of ChlP-PCR, ChlP-Seq with 
bioinformatics analysis, and supplemental figure legends. 

Additional file 4: Table S5. GREAT Gene Ontology terms. This table 
contains the top GO terms assigned by GREAT to clusters 1 to 7 defined 
in Figure 3A. 

Additional file 5: Figure SI. GO term analysis from GREAT for clusters 
1 to 3 and 5 to 7 in Figure 3A. This is related to Figure 3. 

Additional file 6: Figure S8. Validation of Runx2 peaks by ChlP-PCR. 
This figure is related to Figures 4, 5 and 6. 

Additional file 7: Figure S2. Runx2 peaks associated with Bsp gene 
during differentiation. This figure is related to Figure 4. 

Additional file 8: Table S6. Detailed annotation of Runx2 peaks. This 
table is a meta-spreadsheet of Runx2 peaks annotated to RefSeq genes. 

Additional file 9: Figure S3. Validation of Runx2 knockdown in MC3T3 
cells. This figure is related to Figure 5. 

Additional file 10: Figure S4. Additional characteristics of Runx2 
binding in shRunx2 responsive genes. This figure is related to Figure 5. 

Additional file 11: Figure S5. PeaksToGenes analysis of Runx2 occupancy 
in Runx2 shRNA-responsive genes. This figure is related to Figure 5E. 

Additional file 12: Figure S6. EMBER analyses of Runx2 binding in the 
genes differentially regulated by Runx2 knockdown. This figure is related 
to Figure 5. 

Additional file 13: Figure S7. Validation of novel Runx2 target TnfrsfW. 
This figure is related to Figure 6. 

Additional file 14: Table S7. Genes responsive to Runx2 shRNA. This 
table lists the genes that are up- or down- regulated by shRunx2 treat- 
ment in MC3T3 cells differentiated for 9 days. 

Additional file 15: Table SI. qPCR and ChlP-PCR primers. 

Additional file 16: Table S2. Cloning primers. This table contains the 
primers used for plasmid construction. 
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